Mix PCA - Données Géométrie
Madison Giacofci
2023-2024
Chargement des données et preprocessing
- Chargement des données
geom = readxl::read_xlsx("geometriePiece/ResultsAddedForExcel.xlsx",col_names=TRUE,skip=1)New names:
• `` -> `...1`colnames(geom)[1] = "type"
geom$type = as.factor(geom$type)- On se concentre sur un cahoutchouc (le premier P225M909)
geom = split(geom,geom$type)$P225M909
Xgeom = geom[,2:4]
Y = geom[,5:ncol(geom)]
Gap = seq(30,7,by=-0.5)
matplot(Gap,t(Y),type="l",col = "steelblue",lwd=0.25,lty=1,xlab="Gap (mm)",ylab="Efforts (N)")
- on remplace les donnees manquantes par un plateau (correspond à des divergences du modele)
NAcurve = sort(unique(which(is.na(Y),arr.ind = TRUE)[,1]))
for (i in NAcurve){Y[i,which(is.na(Y[i,]))] = Y[i,min(which(is.na(Y[i,])))-1]}
# summary(Y) # OK, plus de donnees manquantes- projection dans une base de splines pour obtenir des coefficients splines (smoothing penalise) sur une base de splines de taille 14 (reduction dimension)
splbasis = create.bspline.basis(rangeval=range(Gap),norder=4,breaks=seq(min(Gap),max(Gap),length=12))
gcv = 1:21
for (i in 1:21){
lambda = exp(i-10)
fdparTemp = fdPar(splbasis,Lfdobj = 2,lambda=lambda)
smoothdata = smooth.basis(argvals=Gap,t(Y),fdParobj = fdparTemp)
gcv[i] = mean(smoothdata$gcv)
}
lambda = exp(which.min(gcv)-10)
fdparTemp = fdPar(splbasis,Lfdobj = 2,lambda=lambda)
smoothdata = smooth.basis(argvals=Gap,t(Y),fdParobj = fdparTemp)
# plotfit.fd(t(Y),Gap,smoothdata$fd,pch=20,cex=0.5) #on est tres proche de l'interpolation
Ycoef = t(smoothdata$fd$coefs)
dim(Ycoef) # 14 coefficients pour n individus[1] 498 14Estimation avec MixPCA
- On estime avec notre modèle mixPCA (ne pas faire tourner)
#data = cbind(Ycoef,Xgeom)
#K=3;maxits=100;tol=1e-4; q = 3;p=14;nx=3
#est0 = estimates(cbind(Ycoef,Xgeom),K=3,par_init=NULL,maxits=300,
# tol=1e-4,
# q = 3,
# p=14,
# nx=3,
# cste = TRUE,
# verbose=TRUE)
#save(est0,file="geometriePiece/resultsGeom-31052023.Rdata")- chargement des résultats (31/05/2023)
load("geometriePiece/resultsGeom-31052023.Rdata")Les résultats sont contenus dans l’objet est0
- Représentation des courbes par groupes
couleur = c("coral1","palegreen3","steelblue")
matplot(Gap,t(Y),type="l",col=couleur[est0$G],lwd=0.4,lty=1,xlab="Gap (mm)",ylab="Efforts (N)")
par(mfrow=c(1,3))
for (i in 1:length(est0$piik)){
matplot(Gap,t(Y[which(est0$G==i),]),type="l",col=couleur[i],ylim=c(0,42),lwd=0.4,lty=1,xlab="Gap (mm)",ylab="Efforts (N)")
}
Canonical correlation analysis
On sélectionne d’abord les individus (courbes) d’un groupe particulier (tel que trouvé par MixPCA). On se concentre ici sur le groupe 1
grpe = 1 # a fixer suivant le groupe sur lequel on se concentre
index.gpe = which(est0$G==grpe)
print(paste("Le groupe ", grpe," contient ", length(index.gpe), " individus"))[1] "Le groupe 1 contient 110 individus"matplot(Gap,t(Y)[,index.gpe],type="l",col="steelblue",lty=1,lwd=0.4,xlab="Gap (mm)",ylab="Efforts (N)")
#reconstruction des courbes moyennes fonctionnelles
splbasis = create.bspline.basis(rangeval=range(Gap),norder=4,breaks=seq(min(Gap),max(Gap),length=12))
mufct = sapply(est0$mu,FUN=function(x) {ft =fd(x,splbasis) ; return(eval.fd(Gap,ft))})
#ajout moyenne groupe 2
lines(Gap,mufct[,grpe],col="red",lwd=1)
On prend les courbes du groupe choisi, que l’on projette sur notre base de splines, ainsi que les géométries associées. On standardise les deux matrices (est-ce qu’il faut centrer uniquement?)
YcoefG = scale(Ycoef[index.gpe,])
XgeomG = scale(Xgeom[index.gpe,])Les données de geometrie contiennent, dans l’ordre : * Thickness * Height * Shape
Analyse des correlations sur les deux groupes de données puis globalement
crosscor = matcor(YcoefG,XgeomG)
img.matcor(crosscor,type=2)
Analyse des correlations canoniques - groupe 1
On se place dans un modèle du type
\[ Y = \mu_{Y} + Q_{Y} \alpha + \varepsilon \]
et
\[ X = \mu_{X} + Q_{X} \alpha + \eta \]
cc1 <- cc(XgeomG,YcoefG)
barplot(cc1$cor, xlab = "Dimension", ylab = "Canonical correlations", ylim = c(0,1))
cc1$cor[1] 0.9987295 0.9793919 0.9360975Les 3 corrélations canoniques sont plutot proches de 1
# Affichage des corrélations canoniques
cor(cc1$scores$xscores,cc1$scores$yscores) [,1] [,2] [,3]
[1,] 9.987295e-01 -4.585389e-16 -5.155542e-16
[2,] 7.673963e-18 9.793919e-01 -6.345879e-17
[3,] 9.007057e-16 9.305821e-15 9.360975e-01dim(cc1$scores$xscores)[1] 110 3Les scores selon X et selon Y (notés \(\alpha\)) sont très fortement corrélés et indépendant
La matrice \(Q_{X}\) est donnée par
cc1$xcoef [,1] [,2] [,3]
Thickness -0.9364973 0.02020858 0.4268935
Height 0.1784473 1.03719849 -0.2303829
Shape 0.2744826 -0.18789032 1.0104856t(cc1$xcoef)%*%cc1$xcoef [,1] [,2] [,3]
[1,] 0.9842113 0.1145874 -0.1635351
[2,] 0.1145874 1.1114919 -0.4201863
[3,] -0.1635351 -0.4201863 1.2563954cc1$xcoef%*%t(cc1$xcoef) Thickness Height Shape
Thickness 1.0596736 -0.2445041 0.1705205
Height -0.2445041 1.1607004 -0.3786974
Shape 0.1705205 -0.3786974 1.1317245xcoef1 = cc1$xcoefCette matrice nous renseigne sur la contribution des variables de géométrie à chaque axe factoriel.
La matrice \(Q_{Y}\) est donnée par
Qy = cc1$ycoef
Qy [,1] [,2] [,3]
bspl4.1 0.04981210 -0.45107308 -0.7668100
bspl4.2 -0.05527291 -0.70640730 -2.3186012
bspl4.3 -0.41249359 -0.27723428 -2.4572998
bspl4.4 -0.26682343 0.99493213 -0.7637618
bspl4.5 -0.14870508 -0.66441240 6.8631132
bspl4.6 -0.15504243 -0.71888441 2.3721017
bspl4.7 0.05636064 1.86620349 -2.4302618
bspl4.8 -0.04081794 0.03082379 -1.4813699
bspl4.9 -0.08066014 -0.19333120 -0.3458932
bspl4.10 -0.06740741 0.19068015 0.2732237
bspl4.11 -0.03838069 0.40307386 0.3060224
bspl4.12 0.13905242 1.31726846 1.2849540
bspl4.13 0.29149592 2.42407301 2.6380758
bspl4.14 0.14690757 1.20261944 1.5097830#t(cc1$ycoef)%*%cc1$ycoef
#V = t(cc1$scores$yscores)%*%YcoefG
#cor(t(V))Elle correspond à nos composantes principales fonctionnelles. En “revenant” dans le domaine fonctionnel, on peut voir si celles-ci sont interprétables en tant que telles.
- Représentation des composantes seules
Qfct = apply(Qy,2,FUN=function(x) {ft =fd(x,splbasis) ; return(eval.fd(Gap,ft))})
Qfct1 = Qfct
matplot(Gap,Qfct,col=c("coral1","palegreen3","steelblue"),type="l")
- Représentation comme deviation autour de la moyenne fonctionnelle (voir Ramsay et Silvermann) :
On représente (à vocation d’interprétation graphique)
\[ \mu_Y(t) \pm \alpha \hat{\psi}_j (t), \qquad j=1,\ldots,q\] où \(\alpha\) peut être choisi plus ou moins “arbitrairement” (normalement, on utilise la valeur propre associée pour moduler le coeff mais ici, ce n’est pas clair pour moi)
for(i in 1:3){
matplot(Gap,t(Y)[,index.gpe],type="l",col="gray",ylab="Efforts",ylim=c(0,42),lwd=0.5)
lines(Gap,mufct[,grpe],col=1)
lines(Gap,mufct[,grpe]+3*Qfct[,i],type="p",pch="+",col=couleur[i])
lines(Gap,mufct[,grpe]-3*Qfct[,i],type="p",pch="-",col=couleur[i],lwd=2)
}
Il n’est pas évident de voir sur cette représentation si les composantes principales contenues dans \(Q_Y\) ont un sens fonctionnel.
- Avec la stratégie varimax
Qvarimax = varimax(Qfct)
par(mfrow=c(1,2))
matplot(Gap,Qfct,col=c("coral1","palegreen3","steelblue"),type="l",main="sans VARIMAX")
matplot(Gap,Qvarimax$loadings,col=c("coral1","palegreen3","steelblue"),type="l",main="avec VARIMAX")
for(i in 1:3){
matplot(Gap,t(Y)[,index.gpe],type="l",col="gray",ylab="Efforts",ylim=c(0,42),lwd=0.5)
lines(Gap,mufct[,grpe],col=1)
lines(Gap,mufct[,grpe]+3*Qvarimax$loadings[,i],type="p",pch="+",col=couleur[i])
lines(Gap,mufct[,grpe]-3*Qvarimax$loadings[,i],type="p",pch="-",col=couleur[i],lwd=2)
}
La différence n’est pas complètement claire…
Par ailleurs, on est limité à 3 composantes principales étant donné qu’on a 3 covariables de géométrie. Philosophiquement, ce n’est sans doute pas étonnant puisqu’on a au final des courbes qui n’ont que 3 degrés de liberté.
Analyse des correlations canoniques - groupe 2
On sélectionne d’abord les individus (courbes) d’un groupe particulier (tel que trouvé par MixPCA). On se concentre ici sur le groupe 2
grpe = 2 # a fixer suivant le groupe sur lequel on se concentre
index.gpe = which(est0$G==grpe)
print(paste("Le groupe ", grpe," contient ", length(index.gpe), " individus"))[1] "Le groupe 2 contient 179 individus"matplot(Gap,t(Y)[,index.gpe],type="l",col="steelblue",lty=1,lwd=0.4,xlab="Gap (mm)",ylab="Efforts (N)")
#reconstruction des courbes moyennes fonctionnelles
splbasis = create.bspline.basis(rangeval=range(Gap),norder=4,breaks=seq(min(Gap),max(Gap),length=12))
mufct = sapply(est0$mu,FUN=function(x) {ft =fd(x,splbasis) ; return(eval.fd(Gap,ft))})
#ajout moyenne groupe 2
lines(Gap,mufct[,grpe],col="red",lwd=1)
On prend les courbes du groupe choisi, que l’on projette sur notre base de splines, ainsi que les géométries associées. On standardise les deux matrices (est-ce qu’il faut centrer uniquement?)
YcoefG = scale(Ycoef[index.gpe,])
XgeomG = scale(Xgeom[index.gpe,])Les données de geometrie contiennent, dans l’ordre : * Thickness * Height * Shape
Analyse des correlations sur les deux groupes de données puis globalement
crosscor = matcor(YcoefG[,1:9],XgeomG)
img.matcor(crosscor,type=2)
Analyse des correlations canoniques
On se place dans un modèle du type
\[ Y = \mu_{Y} + Q_{Y} \alpha + \varepsilon \]
et
\[ X = \mu_{X} + Q_{X} \alpha + \eta \]
cc1 <- cc(XgeomG,YcoefG[,1:9]) # on enlève les 5 derniers coefs (singulrite matrice)
barplot(cc1$cor, xlab = "Dimension", ylab = "Canonical correlations", ylim = c(0,1))
cc1$cor[1] 0.9870316 0.8924068 0.4784289Les 3 corrélations canoniques sont moins proches de 1 ici (interpretation à reflechir)
# Affichage des corrélations canoniques
cor(cc1$scores$xscores,cc1$scores$yscores) [,1] [,2] [,3]
[1,] 9.870316e-01 1.377484e-16 -1.320109e-15
[2,] -4.273553e-17 8.924068e-01 1.458884e-16
[3,] -1.461115e-16 1.758456e-17 4.784289e-01dim(cc1$scores$xscores)[1] 179 3Les scores selon X et selon Y (notés \(\alpha\)) sont très fortement corrélés et indépendant
La matrice \(Q_{X}\) est donnée par
cc1$xcoef [,1] [,2] [,3]
Thickness 0.8410331 -0.7386838 0.3001075
Height 0.4434320 0.9675255 -0.3651595
Shape -0.4890760 0.3397731 0.9320679#t(cc1$xcoef)%*%cc1$xcoef # orthogonalite?
#cc1$xcoef%*%t(cc1$xcoef) # orthogonalite?
xcoef2 = cc1$xcoefCette matrice nous renseigne sur la contribution des variables de géométrie à chaque axe factoriel. Elle semble a priori assez différente de celle du groupe 1 meme si on retrouve un deuxième “axe” lié à la hauteur et un troisième “axe” lié à shape. Note : elle n’est toujours pas orthogonale
La matrice \(Q_{Y}\) est donnée par
Qy = cc1$ycoef
Qy [,1] [,2] [,3]
bspl4.1 0.08938044 -0.52163186 1.4863663
bspl4.2 0.06420464 -0.21975696 -1.7278458
bspl4.3 0.06382005 0.45138160 0.8170389
bspl4.4 0.13458351 -2.38261643 2.3933461
bspl4.5 0.79450371 0.52431956 -2.9886906
bspl4.6 -0.12820481 2.11949191 -0.8735578
bspl4.7 -0.07424761 0.21352924 0.8470028
bspl4.8 0.08737024 -0.08318324 0.3348507
bspl4.9 0.05717563 -0.25571302 -0.1567053#t(cc1$ycoef)%*%cc1$ycoef
#V = t(cc1$scores$yscores)%*%YcoefG
#cor(t(V))Elle correspond à nos composantes principales fonctionnelles. En “revenant” dans le domaine fonctionnel, on peut voir si celles-ci sont interprétables en tant que telles.
- Représentation des composantes seules
Ycoeffin = matrix(apply(YcoefG,2,mean)[10:14],nrow=5,ncol=3) # on reprend la moyenne des coefs pour la fin des courbes (grands gaps)
Qfct = apply(rbind(Qy,Ycoeffin),2,FUN=function(x) {ft =fd(x,splbasis) ; return(eval.fd(Gap,ft))})
Qfct2 = Qfct
matplot(Gap,Qfct,col=c("coral1","palegreen3","steelblue"),type="l")
- Représentation comme deviation autour de la moyenne fonctionnelle (voir Ramsay et Silvermann) :
On représente (à vocation d’interprétation graphique)
\[ \mu_Y(t) \pm \alpha \hat{\psi}_j (t), \qquad j=1,\ldots,q\] où \(\alpha\) peut être choisi plus ou moins “arbitrairement” (normalement, on utilise la valeur propre associée pour moduler le coeff mais ici, ce n’est pas clair pour moi)
for(i in 1:3){
matplot(Gap,t(Y)[,index.gpe],type="l",col="gray",ylab="Efforts",ylim=c(0,42),lwd=0.5)
lines(Gap,mufct[,grpe],col=1)
lines(Gap,mufct[,grpe]+3*Qfct[,i],type="p",pch="+",col=couleur[i])
lines(Gap,mufct[,grpe]-3*Qfct[,i],type="p",pch="-",col=couleur[i],lwd=2)
}
On peut interpréter les courbes en regard de la matrice
xcoef. Ces courbes nous renseigne sur les parties impactées
par la géométrie en fonction de l’axe factoriel considéré. Par exemple,
le deuxième axe de xcoef est positivement lié à
Height mais négativement à Thickness. La
reconstruction fonctionnelle du deuxième axe ycoef est
cohérente avec cela car on observe que les fortes valeurs positives
(grande hauteur/faible épaisseur) montre une croissance de la courbe
plus précoce que pour les fortes valeurs négatives (petites
hauteurs/forte épaisseur)
Analyse des correlations canoniques - groupe 3
On sélectionne d’abord les individus (courbes) d’un groupe particulier (tel que trouvé par MixPCA). On se concentre ici sur le groupe 3
grpe = 3 # a fixer suivant le groupe sur lequel on se concentre
index.gpe = which(est0$G==grpe)
print(paste("Le groupe ", grpe," contient ", length(index.gpe), " individus"))[1] "Le groupe 3 contient 209 individus"matplot(Gap,t(Y)[,index.gpe],type="l",col="steelblue",lty=1,lwd=0.4,xlab="Gap (mm)",ylab="Efforts (N)")
#reconstruction des courbes moyennes fonctionnelles
splbasis = create.bspline.basis(rangeval=range(Gap),norder=4,breaks=seq(min(Gap),max(Gap),length=12))
mufct = sapply(est0$mu,FUN=function(x) {ft =fd(x,splbasis) ; return(eval.fd(Gap,ft))})
#ajout moyenne groupe 2
lines(Gap,mufct[,grpe],col="red",lwd=1)
On prend les courbes du groupe choisi, que l’on projette sur notre base de splines, ainsi que les géométries associées. On voit ici que le groupe 3 n’est pas très homogène, il contient des courbes avec des comportements très différents. De même que pour les deux autres groupes, on standardise les deux matrices (est-ce qu’il faut centrer uniquement?)
YcoefG = scale(Ycoef[index.gpe,])
XgeomG = scale(Xgeom[index.gpe,])Les données de geometrie contiennent, dans l’ordre : * Thickness * Height * Shape
Analyse des correlations sur les deux groupes de données puis globalement
crosscor = matcor(YcoefG,XgeomG)
img.matcor(crosscor,type=2)
Les corrélations croisées présentent un comportement relativement différents dans ce groupe (en 3 blocs au lieu de 2), faisant echo à l’hétérogénéité des courbes dans le groupe.
Analyse des correlations canoniques
On se place dans un modèle du type
\[ Y = \mu_{Y} + Q_{Y} \alpha + \varepsilon \]
et
\[ X = \mu_{X} + Q_{X} \alpha + \eta \]
cc1 <- cc(XgeomG,YcoefG)
barplot(cc1$cor, xlab = "Dimension", ylab = "Canonical correlations", ylim = c(0,1))
cc1$cor[1] 0.9760231 0.9679074 0.3567894Pas de souci de singularité dans ce groupe. La 3eme corrélation canonique est relativement faible ici (interpretation à reflechir)
# Affichage des corrélations canoniques
cor(cc1$scores$xscores,cc1$scores$yscores) [,1] [,2] [,3]
[1,] 9.760231e-01 -1.909368e-16 -5.523116e-16
[2,] 6.306779e-18 9.679074e-01 4.844452e-16
[3,] 1.436371e-16 3.311742e-17 3.567894e-01#dim(cc1$scores$xscores)Les scores selon X et selon Y (notés \(\alpha\)) sont très fortement corrélés et indépendant (sauf la 3eme)
La matrice \(Q_{X}\) est donnée par
cc1$xcoef [,1] [,2] [,3]
Thickness 0.6909726 0.6703926 0.4054039
Height -0.4893503 0.9174937 0.0650506
Shape -0.2755391 -0.1476526 0.9551033#t(cc1$xcoef)%*%cc1$xcoef # orthogonalite?
#cc1$xcoef%*%t(cc1$xcoef) # orthogonalite?
xcoef3 = cc1$xcoefCette matrice nous renseigne sur la contribution des variables de géométrie à chaque axe factoriel. On retrouve une matrice légèrement différente des autres groupes même si l’axe 2 est toujours lié positivement à la hauteur, l’axe 1 est lié à l’épaisseur positivement et négativement à la hauteur. L’axe 3 reste lié à shape.
La matrice \(Q_{Y}\) est donnée par
Qy = cc1$ycoef
Qy [,1] [,2] [,3]
bspl4.1 0.94397278 -0.150832261 -1.147420132
bspl4.2 -0.42956322 -0.807084948 -0.783172061
bspl4.3 0.22854101 0.414421492 1.654343278
bspl4.4 0.57389907 0.753557408 0.498022864
bspl4.5 -0.46776895 0.157110102 0.009861507
bspl4.6 -0.51906239 0.042890446 2.158460635
bspl4.7 -0.26390244 0.033380575 -0.464832741
bspl4.8 0.63205519 0.269904302 -4.024685529
bspl4.9 -0.05048885 0.284290003 1.711892056
bspl4.10 -0.10579046 0.008809529 0.471992567
bspl4.11 -0.03258642 0.109354682 0.313024771
bspl4.12 -0.01188870 -0.086751591 2.634823744
bspl4.13 0.03706895 -0.210887782 3.993030181
bspl4.14 0.04946376 -0.093288446 1.255825534#t(cc1$ycoef)%*%cc1$ycoef
#V = t(cc1$scores$yscores)%*%YcoefG
#cor(t(V))Elle correspond à nos composantes principales fonctionnelles. En “revenant” dans le domaine fonctionnel, on peut voir si celles-ci sont interprétables en tant que telles.
- Représentation des composantes seules
Qfct = apply(Qy,2,FUN=function(x) {ft =fd(x,splbasis) ; return(eval.fd(Gap,ft))})
Qfct3 = Qfct
matplot(Gap,Qfct,col=c("coral1","palegreen3","steelblue"),type="l")
- Représentation comme deviation autour de la moyenne fonctionnelle (voir Ramsay et Silvermann) :
On représente (à vocation d’interprétation graphique)
\[ \mu_Y(t) \pm \alpha \hat{\psi}_j (t), \qquad j=1,\ldots,q\] où \(\alpha\) peut être choisi plus ou moins “arbitrairement” (normalement, on utilise la valeur propre associée pour moduler le coeff mais ici, ce n’est pas clair pour moi)
for(i in 1:3){
matplot(Gap,t(Y)[,index.gpe],type="l",col="gray",ylab="Efforts",ylim=c(0,20),lwd=0.5)
lines(Gap,mufct[,grpe],col=1)
lines(Gap,mufct[,grpe]+3*Qfct[,i],type="p",pch="+",col=couleur[i])
lines(Gap,mufct[,grpe]-3*Qfct[,i],type="p",pch="-",col=couleur[i],lwd=2)
}
Le même type d’interprétation peut prévaloir ici mais les résultats sont plus difficiles à interpréter, le groupe 3 semblant rassembler des courbes très diverses (groupe un peu fourre-tout)
Comparaison entre les 3 groupes
par(mfrow=c(1,3))
image.plot(matrix(rev(as.vector(t(xcoef1))),3,3)[3:1,], main="Groupe 1",xaxt="n",yaxt="n")
axis(side=2,at=c(0,0.5,1),labels=c("Shape","Height","Thickness"),tick=F)
axis(side=1,at=c(0,0.5,1),labels=c("Axe1","Axe2","Axe3"),tick=F)
image.plot(matrix(rev(as.vector(t(xcoef2))),3,3)[3:1,], main="Groupe 2",xaxt="n",yaxt="n")
axis(side=2,at=c(0,0.5,1),labels=c("Shape","Height","Thickness"),tick=F)
axis(side=1,at=c(0,0.5,1),labels=c("Axe1","Axe2","Axe3"),tick=F)
image.plot(matrix(rev(as.vector(t(xcoef3))),3,3)[3:1,], main="Groupe 3",xaxt="n",yaxt="n")
axis(side=2,at=c(0,0.5,1),labels=c("Shape","Height","Thickness"),tick=F)
axis(side=1,at=c(0,0.5,1),labels=c("Axe1","Axe2","Axe3"),tick=F)
On a des matrices relativement différentes dans les 3 groupes, même si on les regarde à transposition près. NB : dans le groupe 2, la CCA est appliquée en enlevant les 5 derniers coefficients splines des courbes.
Au niveau des axes principaux “fonctionnels”, on a aussi des axes assez différents
par(mfrow=c(1,3))
matplot(Gap,Qfct1,col=c("coral1","palegreen3","steelblue"),type="l")
matplot(Gap,Qfct2,col=c("coral1","palegreen3","steelblue"),type="l")
matplot(Gap,Qfct3,col=c("coral1","palegreen3","steelblue"),type="l")
Il est plus difficile d’expliquer les différences entre les groupes dans cette analyse double. Mais ici, on est parti avec des groupes définis à l’avance, il faudra voir ce que cela donne si on construit des groupes en se basant sur la CCA (ou proche)